section6
zhao hongqiang
6/1/2021
1. Show that a \(3 \times 5\) MA is equivalent to a 7 -term weighted moving average with weights of \(0.067,0.133,0.200,0.200,0.200,0.133\), and \(0.067 .\)
3X5 MA = 0.067,0.133,0.200,0.200,0.200,0.133,0.067 WMA
2. The plastics data set consists of the monthly sales (in thousands) of product A for a plastics manufacturer for five years.
- a:Plot the time series of sales of product A. Can you identify seasonal fluctuations and/or a trend-cycle?
- b:Use a classical multiplicative decomposition to calculate the trend-cycle and seasonal indices.
- c:Do the results support the graphical interpretation from part a?
- d:Compute and plot the seasonally adjusted data.
- e:Change one observation to be an outlier (e.g., add 500 to one observation), and recompute the seasonally adjusted data. What is the effect of the outlier?
- d:Does it make any difference if the outlier is near the end rather than in the middle of the time series?
# a.
autoplot(plastics)
# The plot of the data shows that there are seasonal fluctuations and upward trend.
# b.
decompose_plastics <- decompose(plastics,
type = "multiplicative")
autoplot(decompose_plastics)
# c.
# Yes.
# d.
autoplot(plastics, series="Data") +
autolayer(trendcycle(decompose_plastics), series="Trend") +
autolayer(seasadj(decompose_plastics), series="Seasonally Adjusted") +
xlab("Year") + ylab("Monthly Sales amount") +
ggtitle("Sales of plastic product (in thousand)") +
scale_colour_manual(values=c("gray","blue","red"),
breaks=c("Data","Seasonally Adjusted","Trend"))## Warning: Removed 12 row(s) containing missing values (geom_path).
# e.
plastics_new <- plastics
plastics_new[20] <- plastics_new[20] + 500
decompose_plastics_new <- decompose(
plastics_new,
type = "multiplicative"
)
autoplot(plastics_new, series = "Data") +
autolayer(trendcycle(decompose_plastics_new),
series = "Trend") +
autolayer(seasadj(decompose_plastics_new),
series = "Seasonally Adjusted") +
xlab("Year") +
ylab("Monthly Sales amount") +
ggtitle("Sales of plastic projuct with outlier") +
scale_color_manual(values=c("gray", "blue", "red"),
breaks=c("Data", "Seasonally Adjusted", "Trend"))## Warning: Removed 12 row(s) containing missing values (geom_path).
# The outlier affects trend little, but affects the seasonally adjusted data severely. Seasonally adjusted data have errors like the original data have.
# f.
plastics_new[55] <- plastics_new[55] + 500
decompose_plastics_new <- decompose(
plastics_new,
type = "multiplicative"
)
autoplot(plastics_new, series = "Data") +
autolayer(trendcycle(decompose_plastics_new),
series = "Trend") +
autolayer(seasadj(decompose_plastics_new),
series = "Seasonally Adjusted") +
xlab("Year") + ylab("Monthly Sales amount") +
ggtitle("Sales of plastic projuct with outliers") +
scale_color_manual(values=c("gray", "blue", "red"),
breaks=c("Data", "Seasonally Adjusted", "Trend"))## Warning: Removed 12 row(s) containing missing values (geom_path).
# If an outlier is near the end, the effect to trend decreases.3. Recall your retail time series data (from Exercise 3 in Section 2.10). Decompose the series using X11. Does it reveal any outliers, or unusual features that you had not noticed previously?
# Load retail.xlsx data
library(xlsx)
retail <- read.xlsx("retail.xlsx",
sheetIndex = 1,
startRow = 2)
head(retail)## Series.ID A3349335T A3349627V A3349338X A3349398A A3349468W A3349336V
## 1 1982-04-01 303.1 41.7 63.9 408.7 65.8 91.8
## 2 1982-05-01 297.8 43.1 64.0 404.9 65.8 102.6
## 3 1982-06-01 298.0 40.3 62.7 401.0 62.3 105.0
## 4 1982-07-01 307.9 40.9 65.6 414.4 68.2 106.0
## 5 1982-08-01 299.2 42.1 62.6 403.8 66.0 96.9
## 6 1982-09-01 305.4 42.0 64.4 411.8 62.3 97.5
## A3349337W A3349397X A3349399C A3349874C A3349871W A3349790V A3349556W
## 1 53.6 211.3 94.0 32.7 126.7 178.3 50.4
## 2 55.4 223.8 105.7 35.6 141.3 202.8 49.9
## 3 48.4 215.7 95.1 32.5 127.6 176.3 48.0
## 4 52.1 226.3 95.3 33.5 128.8 172.6 48.6
## 5 54.2 217.1 82.8 29.4 112.3 169.6 51.3
## 6 53.6 213.4 89.4 32.2 121.6 181.4 49.6
## A3349791W A3349401C A3349873A A3349872X A3349709X A3349792X A3349789K
## 1 22.2 43.0 62.4 178.0 61.8 85.4 147.2
## 2 23.1 45.3 63.1 181.5 60.8 84.8 145.6
## 3 22.8 43.7 59.6 174.1 58.7 80.7 139.4
## 4 23.2 46.5 61.9 180.2 60.3 82.4 142.7
## 5 21.4 44.8 60.7 178.1 56.1 80.7 136.8
## 6 21.8 43.9 61.2 176.5 58.1 82.1 140.2
## A3349555V A3349565X A3349414R A3349799R A3349642T A3349413L A3349564W
## 1 1250.2 257.9 17.3 34.9 310.2 58.2 55.8
## 2 1300.0 257.4 18.1 34.6 310.1 62.0 58.4
## 3 1234.2 261.2 18.1 34.6 313.9 53.8 53.7
## 4 1265.0 266.1 18.9 35.2 320.2 57.9 56.9
## 5 1217.6 247.2 19.0 33.8 300.1 59.2 56.7
## 6 1244.9 262.4 18.4 35.4 316.2 57.1 58.9
## A3349416V A3349643V A3349483V A3349722T A3349727C A3349641R A3349639C
## 1 59.1 173.1 93.6 26.3 119.9 104.2 42.2
## 2 59.2 179.5 95.3 27.1 122.5 110.2 42.1
## 3 59.8 167.3 85.2 24.3 109.6 96.7 38.5
## 4 59.8 174.5 91.6 25.6 117.2 104.6 38.9
## 5 62.2 178.1 85.2 23.5 108.7 92.5 39.5
## 6 63.6 179.6 89.5 24.3 113.8 98.3 41.7
## A3349415T A3349349F A3349563V A3349350R A3349640L A3349566A A3349417W
## 1 15.6 31.6 34.4 123.7 36.4 48.7 85.1
## 2 15.8 31.5 34.4 123.9 36.2 48.9 85.1
## 3 15.2 29.6 33.5 116.8 35.7 47.1 82.8
## 4 15.2 35.2 33.4 122.7 34.6 47.5 82.1
## 5 14.5 34.7 33.2 122.0 32.5 49.3 81.8
## 6 15.1 34.2 34.5 125.5 33.9 50.7 84.6
## A3349352V A3349882C A3349561R A3349883F A3349721R A3349478A A3349637X
## 1 916.2 139.3 NA NA 161.8 31.8 46.6
## 2 931.2 136.0 NA NA 158.7 32.8 49.6
## 3 887.0 143.5 NA NA 166.6 34.9 51.4
## 4 921.3 150.2 NA NA 172.9 34.6 50.9
## 5 883.2 144.0 NA NA 165.9 32.9 51.6
## 6 917.9 146.9 NA NA 169.5 33.7 49.6
## A3349479C A3349797K A3349477X A3349719C A3349884J A3349562T A3349348C
## 1 13.3 91.6 28.9 13.9 42.8 67.5 18.4
## 2 12.7 95.0 30.6 14.7 45.3 69.7 17.7
## 3 12.9 99.2 30.5 14.5 45.1 60.7 17.7
## 4 13.9 99.4 27.9 15.2 43.1 67.9 18.4
## 5 12.8 97.3 27.4 14.1 41.5 66.5 17.8
## 6 14.5 97.9 29.1 15.5 44.5 73.4 18.8
## A3349480L A3349476W A3349881A A3349410F A3349481R A3349718A A3349411J
## 1 11.1 22.0 25.8 77.3 18.7 26.7 45.4
## 2 11.7 21.9 25.9 77.2 19.5 27.3 46.8
## 3 11.5 22.7 25.9 77.7 18.6 26.2 44.8
## 4 13.1 24.3 28.7 84.4 22.6 25.2 47.8
## 5 13.0 23.6 27.7 82.1 22.6 25.6 48.2
## 6 13.0 21.8 29.0 82.6 23.2 26.7 49.8
## A3349638A A3349654A A3349499L A3349902A A3349432V A3349656F A3349361W
## 1 486.3 83.5 6.0 11.3 100.8 15.2 16.0
## 2 492.8 80.6 5.4 11.1 97.1 17.2 19.0
## 3 494.1 82.3 5.2 11.2 98.7 17.4 18.1
## 4 515.6 88.2 5.6 12.1 105.9 18.7 20.3
## 5 501.4 82.3 5.7 11.7 99.7 18.6 19.6
## 6 517.7 84.2 5.8 12.0 102.0 18.8 19.9
## A3349501L A3349503T A3349360V A3349903C A3349905J A3349658K A3349575C
## 1 8.6 39.7 19.1 6.6 25.7 48.9 8.1
## 2 9.5 45.7 21.6 7.0 28.6 52.2 7.5
## 3 8.4 43.9 18.3 6.0 24.3 48.9 6.7
## 4 10.3 49.3 18.6 6.4 25.0 48.3 7.8
## 5 10.6 48.9 17.1 6.0 23.1 49.4 7.9
## 6 11.5 50.2 18.2 6.4 24.6 48.5 7.8
## A3349428C A3349500K A3349577J A3349433W A3349576F A3349574A A3349816F
## 1 6.1 7.2 12.9 34.2 14.3 15.8 30.1
## 2 6.5 7.5 13.0 34.4 14.2 15.8 30.0
## 3 6.1 7.5 12.5 32.7 13.4 15.3 28.7
## 4 6.6 7.9 13.9 36.2 14.5 17.0 31.4
## 5 6.3 8.3 13.7 36.1 13.6 17.5 31.1
## 6 6.4 7.8 14.1 36.0 13.9 17.8 31.7
## A3349815C A3349744F A3349823C A3349508C A3349742A A3349661X A3349660W
## 1 279.4 96.6 12.3 13.1 122.0 19.2 22.5
## 2 288.0 96.4 11.8 13.4 121.6 21.9 27.8
## 3 277.2 95.6 11.3 13.5 120.4 19.9 26.7
## 4 296.1 103.3 12.1 13.8 129.2 19.3 28.2
## 5 288.4 96.6 12.0 13.3 121.9 19.6 27.4
## 6 293.0 101.4 12.3 13.4 127.1 19.9 27.0
## A3349909T A3349824F A3349507A A3349580W A3349825J A3349434X A3349822A
## 1 8.6 50.4 21.4 7.4 28.8 36.5 9.7
## 2 8.2 57.9 24.1 8.0 32.1 43.7 11.0
## 3 7.9 54.4 21.4 7.0 28.5 38.0 10.7
## 4 8.7 56.2 21.8 7.2 29.0 42.0 9.0
## 5 7.9 55.0 18.7 6.6 25.3 38.5 9.1
## 6 8.7 55.6 19.5 7.4 26.9 40.2 10.0
## A3349821X A3349581X A3349908R A3349743C A3349910A A3349435A A3349365F
## 1 6.5 14.6 11.3 42.1 8.0 10.4 18.4
## 2 7.2 15.2 11.6 45.0 8.0 10.3 18.3
## 3 6.6 14.5 10.9 42.5 7.3 10.4 17.7
## 4 7.0 14.6 11.4 42.0 7.8 10.3 18.1
## 5 6.8 15.3 10.9 42.1 7.6 10.1 17.7
## 6 7.1 15.1 11.7 43.9 8.2 10.3 18.5
## A3349746K A3349370X A3349754K A3349670A A3349764R A3349916R A3349589T
## 1 298.3 26.0 NA NA 28.4 6.1 5.1
## 2 318.5 25.4 NA NA 27.7 6.3 4.7
## 3 301.5 25.3 NA NA 27.7 6.4 5.2
## 4 316.4 27.8 NA NA 30.3 5.9 5.2
## 5 300.5 26.6 NA NA 29.0 5.7 4.8
## 6 312.3 27.1 NA NA 29.6 5.3 4.8
## A3349590A A3349765T A3349371A A3349588R A3349763L A3349372C A3349442X
## 1 2.4 13.6 6.7 1.9 8.7 NA 2.9
## 2 2.5 13.4 7.4 1.9 9.3 NA 2.9
## 3 2.1 13.7 6.7 1.8 8.6 NA 2.9
## 4 2.7 13.7 7.1 1.8 8.9 NA 3.1
## 5 2.9 13.4 5.8 1.7 7.5 NA 3.1
## 6 2.6 12.8 5.8 1.7 7.5 NA 3.2
## A3349591C A3349671C A3349669T A3349521W A3349443A A3349835L A3349520V
## 1 1.8 4.0 NA NA 1.9 3.5 5.4
## 2 1.9 4.0 NA NA 2.0 3.5 5.5
## 3 1.9 3.9 NA NA 2.0 3.1 5.1
## 4 1.8 4.4 NA NA 1.9 3.6 5.5
## 5 1.8 4.2 NA NA 1.9 3.6 5.5
## 6 1.8 4.0 NA NA 1.9 3.8 5.7
## A3349841J A3349925T A3349450X A3349679W A3349527K A3349526J A3349598V
## 1 79.9 NA NA NA NA NA NA
## 2 78.9 NA NA NA NA NA NA
## 3 77.5 NA NA NA NA NA NA
## 4 82.7 NA NA NA NA NA NA
## 5 78.1 NA NA NA NA NA NA
## 6 79.1 NA NA NA NA NA NA
## A3349766V A3349600V A3349680F A3349378T A3349767W A3349451A A3349924R
## 1 NA NA NA NA NA NA NA
## 2 NA NA NA NA NA NA NA
## 3 NA NA NA NA NA NA NA
## 4 NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA
## A3349843L A3349844R A3349376L A3349599W A3349377R A3349779F A3349379V
## 1 NA NA NA NA NA NA NA
## 2 NA NA NA NA NA NA NA
## 3 NA NA NA NA NA NA NA
## 4 NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA
## A3349842K A3349532C A3349931L A3349605F A3349688X A3349456L A3349774V
## 1 NA 12.7 1.2 1.6 15.5 2.7 4.4
## 2 NA 12.1 1.4 1.6 15.1 3.0 4.9
## 3 NA 12.5 1.3 1.7 15.5 2.5 4.8
## 4 NA 13.2 1.4 1.6 16.1 2.8 5.1
## 5 NA 12.7 1.6 1.6 15.8 2.8 4.6
## 6 NA 12.9 1.4 1.8 16.0 2.6 4.3
## A3349848X A3349457R A3349851L A3349604C A3349608L A3349609R A3349773T
## 1 2.6 9.7 3.7 2.2 5.9 10.3 2.3
## 2 3.3 11.1 3.8 2.1 5.9 10.6 2.5
## 3 2.7 9.9 3.2 2.0 5.1 9.9 2.3
## 4 2.4 10.2 3.4 2.1 5.4 8.8 2.6
## 5 2.7 10.1 3.1 2.0 5.0 8.8 2.6
## 6 3.1 10.0 3.4 2.2 5.6 9.2 2.6
## A3349852R A3349775W A3349776X A3349607K A3349849A A3349850K A3349606J
## 1 1.1 2.5 2.2 8.1 4.4 3.2 7.6
## 2 1.0 2.5 2.0 8.0 3.4 3.3 6.7
## 3 1.0 2.5 2.0 7.8 3.6 3.5 7.1
## 4 1.1 2.6 2.0 8.3 4.0 3.5 7.5
## 5 0.9 2.8 2.0 8.4 3.6 3.7 7.3
## 6 1.0 2.8 2.2 8.6 4.2 3.9 8.1
## A3349932R A3349862V A3349462J A3349463K A3349334R A3349863W A3349781T
## 1 57.1 933.4 79.6 149.6 1162.6 200.3 243.4
## 2 57.3 920.5 80.8 149.7 1150.9 210.3 268.3
## 3 55.3 933.6 77.3 149.0 1160.0 198.7 266.1
## 4 56.3 972.6 80.4 153.5 1206.4 208.7 273.5
## 5 55.4 923.5 81.6 147.3 1152.5 206.2 262.7
## 6 57.5 955.9 81.4 151.8 1189.1 200.9 263.1
## A3349861T A3349626T A3349617R A3349546T A3349787F A3349333L A3349860R
## 1 148.6 592.3 268.5 91.4 359.9 460.1 135.1
## 2 151.0 629.6 289.8 96.8 386.6 502.6 134.9
## 3 142.6 607.4 261.9 88.6 350.5 443.8 128.2
## 4 150.1 632.4 267.2 92.1 359.3 459.1 129.9
## 5 153.7 622.6 241.5 83.7 325.2 438.4 133.0
## 6 157.9 622.0 256.2 90.1 346.3 465.1 135.5
## A3349464L A3349389X A3349461F A3349788J A3349547V A3349388W A3349870V
## 1 64.9 125.6 153.5 479.1 146.3 196.1 342.4
## 2 67.7 128.7 154.8 486.1 145.5 196.6 342.1
## 3 65.5 125.0 148.8 467.5 140.2 188.5 328.7
## 4 68.5 136.6 156.1 491.1 146.5 192.0 338.5
## 5 65.2 134.7 152.8 485.7 138.8 192.7 331.5
## 6 66.8 130.4 157.2 489.9 144.3 197.6 341.9
## A3349396W NA. NA..1 NA..2 NA..3 NA..4 NA..5 NA..6 NA..7 NA..8 NA..9
## 1 3396.4 NA NA NA NA NA NA NA NA NA NA
## 2 3497.9 NA NA NA NA NA NA NA NA NA NA
## 3 3357.8 NA NA NA NA NA NA NA NA NA NA
## 4 3486.8 NA NA NA NA NA NA NA NA NA NA
## 5 3355.9 NA NA NA NA NA NA NA NA NA NA
## 6 3454.3 NA NA NA NA NA NA NA NA NA NA# change the retail data to ts object
ts_retail <- ts(retail[,"A3349873A"],
frequency=12,
start=c(1982,4))
# plot the data to see if there is a trend or seasonality.
autoplot(ts_retail)
# There are trend and seasonality.
# Decompose the series using X11.
x11_retail <- seas(ts_retail, x11 = "")
autoplot(x11_retail)
# there were some outliers that I didn't expect. Especially the biggest outlier happened in 2001. And I didn't expect the seasonality effect decreases as trend increases.4. The result of decomposing the number of persons in the civilian labor force in Australia each month from February 1978 to August 1995 are shown below.
# First figure: Decomposition of the data of the number of persons in the civilian labor force in Australia each month from February 1978 to August 1995.
library(seasonal)
library(fpp2)
stl_labour <- stl(labour, s.window = 13, robust = TRUE)
autoplot(stl_labour) + xlab("Year") +
ggtitle("STL decomposition of civilian labour force")
# Second figure: ggsubseries plot of civilian labor force.
ggsubseriesplot(seasonal(stl_labour)) +
ylab("Seasonal")
- a:Write about 3–5 sentences describing the results of the decomposition. Pay particular attention to the scales of the graphs in making your interpretation.
- b:Is the recession of 1991/1992 visible in the estimated components?
# a.
autoplot(labour, series="Data") +
autolayer(trendcycle(stl_labour), series="Trend") +
autolayer(seasadj(stl_labour), series="Seasonally Adjusted") +
xlab("Year") + ylab("Number of people") +
ggtitle("Number of people in the civilian labour force in Australia") +
scale_colour_manual(values=c("gray","blue","red"),
breaks=c("Data","Seasonally Adjusted","Trend"))
# Overall, the number of people in the civilian labour force in Australia increased over time. There were big recessions around 1991 and 1992, and seasonally adjusted data reflected the events. From ggsubseriesplot, I could find that there were relatively small changes for seasonal component over time relative to the original data scale.
# b.
# It is seen in the seasonally adjuested data.- This exercise uses the cangas data (monthly Canadian gas production in billions of cubic metres, January 1960 - February 2005).
- a: Plot the data using autoplot(), ggsubseriesplot() and ggseasonplot() to look at the effect of the changing seasonality over time. What do you think is causing it to change so much?
- b:Do an STL decomposition of the data. You will need to choose s.window to allow for the changing shape of the seasonal component.
- c:Compare the results with those obtained using SEATS and X11. How are they different?
# a.
autoplot(cangas)
ggsubseriesplot(cangas)
ggseasonplot(cangas)
# The gas production amount increased in winter and decreased in summer. Maybe the cold weather in winter increased the demand for the gas that the production amount increased. But the production amount also increased in summer as time went on. Maybe it happended because the demand for electricity to run air conditioners increased over time.
# b.
stl_cangas <- stl(cangas, s.window = 13, robust = TRUE)
# Show each STL decomposed component
autoplot(stl_cangas) +
ggtitle("Monthly Canadian Gas Production",
subtitle = "STL decomposition")
# can see that the size of seasonality increased in 1980s and decreased in 1990s.
# (STL decomposition) See trend-cycle component and seasonally adjusted data along with original data.
autoplot(cangas, series = "Data") +
autolayer(seasadj(stl_cangas), series = "Seasonally Adjusted") +
autolayer(trendcycle(stl_cangas), series = "Trend-cycle") +
ggtitle("Monthly Canadian gas production(STL decomposition)") +
ylab(expression(paste("Gas production (x", 10^{9}, m^{3}, ")"))) +
scale_color_manual(values = c("gray", "blue", "red"),
breaks = c("Data", "Seasonally Adjusted", "Trend-cycle"))
# c.
x11_cangas <- seas(cangas, x11 = "")
seats_cangas <- seas(cangas)
# Show each X11 decomposed component
autoplot(x11_cangas) +
ggtitle("Monthly Canadian Gas Production",
subtitle = "X11 decomposition")
# (X11 decomposition) See trend-cycle component and seasonally adjusted data along with original data.
autoplot(cangas, series = "Data") +
autolayer(seasadj(x11_cangas), series = "Seasonally Adjusted") +
autolayer(trendcycle(x11_cangas), series = "Trend-cycle") +
ggtitle("Monthly Canadian gas production(X11 decomposition)") +
ylab(expression(paste("Gas production (x", 10^{9}, m^{3}, ")"))) +
scale_color_manual(values = c("gray", "blue", "red"),
breaks = c("Data", "Seasonally Adjusted", "Trend-cycle"))
# Show each SEATS decomposed component
autoplot(seats_cangas) +
ggtitle("Monthly Canadian Gas Production",
subtitle = "SEATS decomposition")
# (SEATS decomposition) See trend-cycle component and seasonally adjusted data along with original data.
autoplot(cangas, series = "Data") +
autolayer(seasadj(seats_cangas), series = "Seasonally Adjusted") +
autolayer(trendcycle(seats_cangas), series = "Trend-cycle") +
ggtitle("Monthly Canadian gas production(SEATS decomposition)") +
ylab(expression(paste("Gas production (x", 10^{9}, m^{3}, ")"))) +
scale_color_manual(values = c("gray", "blue", "red"),
breaks = c("Data", "Seasonally Adjusted", "Trend-cycle"))
# seas function did multiplicative decomposition. Therefore seasonal component and remainder component have mean at 1, not 0. And the proportion of seasonality to trend decreased, then increased, and then decreased again.
# From the plots, I could see that the seasonally adjusted data of STL decomposition have higher variance than the other methods. It can mean that STL method is more robust, that is, unusual observations affect the remainder component more when STL method is used.- We will use the bricksq data (Australian quarterly clay brick production. 1956-1994) for this exercise.
- a:Use an STL decomposition to calculate the trend-cycle and seasonal indices. (Experiment with having fixed or changing seasonality.)
- b:Compute and plot the seasonally adjusted data.
- c:Use a naïve method to produce forecasts of the seasonally adjusted data.
- d:Use stlf() to reseasonalise the results, giving forecasts for the original data.
- e:Do the residuals look uncorrelated?
- f:Repeat with a robust STL decomposition. Does it make much difference?
- g:Compare forecasts from stlf() with those from snaive(), using a test set comprising the last 2 years of data. Which is better?
# a.
# STL decomposition with fixed seasonality
stl_brick_fixed_st <- stl(bricksq,
s.window = "periodic",
robust = TRUE)
# STL decomposition with changing seasonality
stl_brick_changing_st <- stl(bricksq,
s.window = 5,
robust = TRUE)
# plot decomposed data
autoplot(stl_brick_fixed_st) +
ggtitle("Brick production data decomposed by STL with fixed seasonality")
autoplot(stl_brick_changing_st) +
ggtitle("Brick production data decomposed by STL with changing seasonality")
# can see changing seasonal component and smaller remainders.
# b.
# plot data which are decomposed by STL with fixed seasonality
autoplot(bricksq, series = "Data") +
autolayer(trendcycle(stl_brick_fixed_st),
series = "Trend-cycle") +
autolayer(seasadj(stl_brick_fixed_st),
series = "Seasonally Adjusted Data") +
ggtitle("Quarterly clay brick production in Australia",
subtitle = "-decomposed by STL with fixed seasonality") +
scale_color_manual(values = c("gray", "red", "blue"),
breaks = c("Data", "Trend-cycle", "Seasonally Adjusted Data"))
# plot data which are decomposed by STL with changing seasonality
autoplot(bricksq, series = "Data") +
autolayer(trendcycle(stl_brick_fixed_st),
series = "Trend-cycle") +
autolayer(seasadj(stl_brick_fixed_st),
series = "Seasonally Adjusted Data") +
ggtitle("Quarterly clay brick production in Australia",
subtitle = "-decomposed by STL with changing seasonality") +
scale_color_manual(values = c("gray", "red", "blue"),
breaks = c("Data", "Trend-cycle", "Seasonally Adjusted Data"))
# c.
stl_brick_fixed_st %>% seasadj() %>% naive() %>% autoplot() +
ggtitle(label = "Naive forecast of seasonally adjusted brick data",
subtitle = "after STL decomposition with fixed seasonality")
stl_brick_changing_st %>% seasadj() %>% naive() %>% autoplot() +
ggtitle(label = "Naive forecast of seasonally adjusted brick data",
subtitle = "after STL decomposition with changing seasonality")
# can see that the prediction intervals of seasonally adjusted data decomposed by STL with changing seasonality have smaller range than the one with fixed seasonality. It happened because the variance of the remainder component decreased when the seasonality can be changed.
# d.
stlf_brick <- stlf(bricksq)
autoplot(stlf_brick)
# e.
checkresiduals(stlf_brick)## Warning in checkresiduals(stlf_brick): The fitted degrees of freedom is based on
## the model used for the seasonally adjusted data.
##
## Ljung-Box test
##
## data: Residuals from STL + ETS(M,N,N)
## Q* = 41.128, df = 6, p-value = 2.733e-07
##
## Model df: 2. Total lags used: 8# The residuals are correlated with each other.
# f.
stlf_brick_robust <- stlf(bricksq, robust = TRUE)
autoplot(stlf_brick_robust)
checkresiduals(stlf_brick_robust)## Warning in checkresiduals(stlf_brick_robust): The fitted degrees of freedom is
## based on the model used for the seasonally adjusted data.
##
## Ljung-Box test
##
## data: Residuals from STL + ETS(M,N,N)
## Q* = 28.163, df = 6, p-value = 8.755e-05
##
## Model df: 2. Total lags used: 8# It looked like the autocorrelations became lower generally, but there are still some high values left.
# g.
trainset_brick <- subset(bricksq,
end = length(bricksq) - 8)
testset_brick <- subset(bricksq,
start = length(bricksq) - 7)
snaive_brick <- snaive(trainset_brick)
stlf_brick_part <- stlf(trainset_brick, robust = TRUE)
# plot data and forecast results
autoplot(bricksq, series = "Original data") +
geom_line(size = 1) +
autolayer(stlf_brick_part, PI = FALSE, size = 1,
series = "stlf") +
autolayer(snaive_brick, PI = FALSE, size = 1,
series = "snaive") +
scale_color_manual(values = c("gray50", "blue", "red"),
breaks = c("Original data", "stlf", "snaive")) +
scale_x_continuous(limits = c(1990, 1994.5)) +
scale_y_continuous(limits = c(300, 600)) +
guides(colour = guide_legend(title = "Data")) +
ggtitle("Forecast from stlf and snaive functions") +
annotate(
"rect",
xmin=1992.75,xmax=1994.5,ymin=-Inf,ymax=Inf,
fill="lightgreen",alpha = 0.3
)## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 136 row(s) containing missing values (geom_path).
# can see from the plot that the forecasts from stlf function are more similar to the original data than the forecasts from snaive function. Unlike snaive function, stlf function can also use trend, and its seasonality can change over time. The test set have trend with seasonality. Therefore stlf function was better than snaive function to predict brick production amount of near future.- Use stlf to produce forecasts of the writing series with either method=“naive” or method=“rwdrift”, whichever is most appropriate. Use the lambda argument if you think a Box-Cox transformation is required.
str(writing)## Time-Series [1:120] from 1968 to 1978: 563 599 669 598 580 ...head(writing)## Jan Feb Mar Apr May Jun
## 1968 562.674 599.000 668.516 597.798 579.889 668.233autoplot(writing)
# can see that there are increasing trend in writing data. Therefore it would be better to use rwdrift method to forecast non-seasonal(seasonally adjusted) component.
# I think that it would be better to do Box-Cox transformation to make the size of the seasonal variation in the data about the same across the whole series.
stlf_writing <- stlf(writing,
s.window = 13,
robust = TRUE,
lambda = BoxCox.lambda(writing),
method = "rwdrift")
autoplot(stlf_writing)
- Use stlf to produce forecasts of the fancy series with either method=“naive” or method=“rwdrift”, whichever is most appropriate. Use the lambda argument if you think a Box-Cox transformation is required.
str(fancy)## Time-Series [1:84] from 1987 to 1994: 1665 2398 2841 3547 3753 ...head(fancy)## Jan Feb Mar Apr May Jun
## 1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74autoplot(fancy)
# There is increasing trend that it would be better to use rwdrift method for forecasting non-seasonal component.
# And it would be better to do Box-Cox transformation to make the variation in the data about the same across the whole series.
stlf_fancy <- stlf(fancy,
s.window = 13,
robust = TRUE,
lambda = BoxCox.lambda(fancy),
method = "rwdrift")
autoplot(stlf_fancy)
# The prediction intervals increase dramatically because of Box-Cox transformation. But without the transformation, the forecasts are unreasonable.